///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/fqmc/linear_method_classes.h
///
/// \brief   header file for the linear method optimizer
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_FQMC_LINEAR_METHOD_CLASSES_HEADER
#define FORMIC_FQMC_LINEAR_METHOD_CLASSES_HEADER

#include<vector>
#include<string>
#include<numeric>
#include<cassert>
#include<algorithm>
#include<cmath>
#include<complex>

#include<boost/format.hpp>
#include<boost/shared_ptr.hpp>

#include<formic/fqmc/linear_method.h>
#include<formic/wfn/configuration.h>
#include<formic/fqmc/vmc_walker.h>
#include<formic/fqmc/energy.h>
#include<formic/wfn/wfn.h>
#include<formic/ham/ham.h>
#include<formic/input/base.h>
#include<formic/mpi/interface.h>
#include<formic/numeric/numeric.h>
#include<formic/numeric/matrix.h>
#include<formic/random/random.h>
#include<formic/output/output.h>
#include<formic/exception.h>
#include<formic/lapack/interface.h>
#include<formic/davidson.h>
#include<formic/genminres.h>
#include<formic/vector/vector.h>
#include<formic/vector/matrix.h>
#include<formic/timing/timing.h>
#include<formic/fqmc/sampler.h>

namespace formic { namespace fqmc {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a base class for accumulating the result of a linear transformation in the basis of
  ///         the wavefunction and its first derivatives
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class DerivLTBaseAcc : public formic::Accumulator<S> {

    protected:

//      /// \brief  pointer to the coefficient object to use
//      boost::shared_ptr<formic::CoeffObj<S> > _co;

      /// \brief  the number of derivatives
      int _nderiv;

      /// \brief  pointer to the coefficients of the derivative functions in the sum
      const S * _der_coeff;

      /// \brief  pointer to the vector of projections (into our first der. basis) of the result of the lin. trans.
      S * _vec_proj;

      /// \brief  vector to hold the local derivative ratio vector
      std::vector<S> _der_ratios;

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  computes the local contraction for the linear transformation
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      virtual S local_contraction() = 0;

    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  sets the derivative coefficient vector to use
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void set_dc(const S * der_coeff) { _der_coeff = der_coeff; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  sets the vector in which to place the projections
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void set_vp(S * vec_proj) { _vec_proj = vec_proj; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructs the accumulator around a specific coefficient object
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      DerivLTBaseAcc(const formic::InputBase & userinp,
                     const boost::shared_ptr<formic::CoeffObj<S> > & gf,
                     const boost::shared_ptr<formic::CoeffObj<S> > & vf)
        : formic::Accumulator<S>(gf, vf),
          _nderiv(vf->n_der()),
          _der_ratios(1+vf->n_der(), formic::zero(S()))
      {}

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief Accumulate the local contribution
      ///
      /// \param[in]      userinp  the user's input options
      /// \param[in]      w        the weight to apply to this configuration's contribution
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void child_accumulate(const formic::InputBase & userinp, const double w) {

        // get the local contraction
        const S local_cont = this->local_contraction(); //_co->der_hc(_der_coeff);

        // get the derivative ratios and complex conjugate them
        _der_ratios.at(0) = formic::unity(S());
        this->_base_vf->der_ratios(&_der_ratios.at(0) + 1);
        formic::xconj(_der_ratios.size(), &_der_ratios.at(0), 1);

        // add this configuration's contribution to the vector of projections
        formic::xaxpy(_nderiv+1, w * this->_vgs * local_cont, &_der_ratios.at(0), 1, _vec_proj, 1);

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief finalizes the accumulation (this may include for example communicating the results)
      ///
      /// \param[in]      userinp  the user's input options
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void child_finalize(const formic::InputBase & userinp) {
        formic::xcopy(_nderiv+1, _vec_proj, 1, &_der_ratios.at(0), 1);
        formic::mpi::reduce(&_der_ratios.at(0), _vec_proj, _nderiv+1, MPI::SUM);
        formic::xscal(_nderiv+1, formic::unity(S()) / this->_total_weight / this->_vgsa, _vec_proj, 1);
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  resets the accumulator
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void child_reset(const formic::Configuration & config) {
        formic::xscal(_nderiv+1, formic::zero(S()), _vec_proj, 1);
      }

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class for accumulating the result of multiplication by the Hamiltonian in the basis
  ///         of the wavefunction and its first derivatives
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class DerivMultHamAcc : public formic::fqmc::DerivLTBaseAcc<S> {

    protected:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  computes the local contraction for the linear transformation
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S local_contraction() { return this->_base_vf->der_hc(this->_der_coeff); }

    public:

      /// \brief  return the accumulator's name
      std::string name() const { return "DerivMultHamAcc accumulator"; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructs the accumulator around a specific coefficient object
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      DerivMultHamAcc(const formic::InputBase & userinp,
                      const boost::shared_ptr<formic::CoeffObj<S> > & gf,
                      const boost::shared_ptr<formic::CoeffObj<S> > & vf)
        : DerivLTBaseAcc<S>(userinp, gf, vf)
      {}

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class for accumulating the result of multiplication by the overlap matrix in the
  ///         basis of the wavefunction and its first derivatives
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class DerivMultOvlpAcc : public formic::fqmc::DerivLTBaseAcc<S> {

    protected:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  computes the local contraction for the linear transformation
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      S local_contraction() { return this->_base_vf->der_ic(this->_der_coeff); }

    public:

      /// \brief  return the accumulator's name
      std::string name() const { return "DerivMultOvlpAcc accumulator"; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructs the accumulator around a specific coefficient object
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      DerivMultOvlpAcc(const formic::InputBase & userinp,
                      const boost::shared_ptr<formic::CoeffObj<S> > & gf,
                      const boost::shared_ptr<formic::CoeffObj<S> > & vf)
        : DerivLTBaseAcc<S>(userinp, gf, vf)
      {}

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class to build the hamiltonian and overlap matrices via sampling
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class HamOvlpBuilderLM {

    private:

      /// \brief  the shift to apply to the diagonal of the Hamiltonian
      S _shift;

      /// \brief  the hamiltonian matrix
      formic::Matrix<S,int> _ham;

      /// \brief  the overlap matrix
      formic::Matrix<S,int> _ovl;

      /// \brief  class for an accumulator that will build the matrices during sampling
      class HamOvlpAcc : public formic::Accumulator<S> {

        private:

          /// \brief  weight for variance
          double _vw;

          /// \brief  weight for energy
          double _ew;

          formic::Matrix<S,int> & _href;
          formic::Matrix<S,int> & _sref;
          //boost::shared_ptr<formic::CoeffObj<S> > _gf;

          formic::Vector<S,int> _der_rat;
          formic::Vector<S,int> _le_der;

          /// \brief  cumulative average of the local energy
          S _le_avg;

          /// \brief  cumulative average of the square of the local energy
          S _le2_avg;

          /// \brief  cumulative average of the local energy derivatives
          formic::Vector<S,int> _led_avg;

          /// \brief  cumulative average of the local energy times local energy derivatives
          formic::Vector<S,int> _led_le_avg;

          /// \brief  cumulative average of the derivative ratios
          formic::Vector<S,int> _dr_avg;

          /// \brief  cumulative average of the derivative ratios times the local energy
          formic::Vector<S,int> _dr_le_avg;

          /// \brief  cumulative average of the derivative ratios times the square of the local energy
          formic::Vector<S,int> _dr_le2_avg;

        public:

          HamOvlpAcc(const boost::shared_ptr<formic::CoeffObj<S> > & gf,
                     const boost::shared_ptr<formic::CoeffObj<S> > & vf,
                     formic::Matrix<S,int> & h,
                     formic::Matrix<S,int> & s,
                     const double var_weight)
          : formic::Accumulator<S>(gf, vf),
            _vw(var_weight),
            _ew(1.0 - var_weight),
            _href(h),
            _sref(s),
            _der_rat(1+vf->n_der(), formic::zero(S())),
            _le_der(1+vf->n_der(), formic::zero(S())),
            _le_avg(formic::zero(S())),
            _le2_avg(formic::zero(S())),
            _led_avg(1+vf->n_der(), formic::zero(S())),
            _led_le_avg(1+vf->n_der(), formic::zero(S())),
            _dr_avg(1+vf->n_der(), formic::zero(S())),
            _dr_le_avg(1+vf->n_der(), formic::zero(S())),
            _dr_le2_avg(1+vf->n_der(), formic::zero(S()))
          {}

          /// \brief  return the accumulator's name
          std::string name() const { return "HamOvlpAcc accumulator"; }

          ///////////////////////////////////////////////////////////////////////////////////////////////////
          /// \brief Accumulate the local contributions to the matrices
          ///
          /// \param[in]      userinp  the user's input options
          /// \param[in]      w        the weight to apply to this configuration's contribution
          ///
          ///////////////////////////////////////////////////////////////////////////////////////////////////
          void child_accumulate(const formic::InputBase & userinp, const double w) {

            // get derivative ratios
            _der_rat(0) = formic::zero(S());
            this->_base_vf->der_ratios(&_der_rat.at(1));

            // get local energy and its derivatives
            this->_base_vf->local_energy_derivatives(_le_der(0), &_le_der.at(1));
            const S le = _le_der(0);
            _le_der(0) = formic::zero(S());

            // include the value to guiding square ratio in the weight
            const S ww = w * this->_vgs * formic::unity(S());

            // add contributions to accumulated constants and vectors
            _le_avg += ww * le; // average of local energies
            _le2_avg += ww * le * le; // average of squares of local energies
            if ( _vw > 0.0 ) {
              _led_avg.axpy(ww, _le_der); // average of local energy derivatives
              _led_le_avg.axpy(ww * le, _le_der); // average of local energy derivatives times local energy
              _dr_avg.axpy(ww, _der_rat); // average of derivative ratios
              _dr_le_avg.axpy(ww * le, _der_rat); // average of derivative ratios times local energy
              _dr_le2_avg.axpy(ww * le * le, _der_rat); // average of derivative ratios times square of local energy
            }

            // get matrix dimension
            const int m = _le_der.size();

            // add variance contribution to hamiltonian matrix
            if ( _vw > 0.0 )
              formic::xgemm('N', 'N', m, m, 1, ww * _vw, &_le_der.at(0), m, &_le_der.at(0), 1, formic::unity(S()), &_href(0,0), m);

            // get "derivative local energies" ( i.e. <n|H|Psi^x> / <n|Psi> )
            formic::xaxpy(m-1, le, &_der_rat.at(1), 1, &_le_der.at(1), 1);

            // add energy contribution to hamiltonian matrix
            _der_rat(0) = formic::unity(S());
            _le_der(0) = le;
            formic::xgemm('N', 'N', m, m, 1, ww * _ew, &_der_rat.at(0), m, &_le_der.at(0), 1, formic::unity(S()), &_href(0,0), m);

            // add overlap matrix contribution
            formic::xgemm('N', 'N', m, m, 1, ww, &_der_rat.at(0), m, &_der_rat.at(0), 1, formic::unity(S()), &_sref(0,0), m);

          }

          ///////////////////////////////////////////////////////////////////////////////////////////////////
          /// \brief finalizes the accumulation (this may include for example communicating the results)
          ///
          /// \param[in]      userinp  the user's input options
          ///
          ///////////////////////////////////////////////////////////////////////////////////////////////////
          void child_finalize(const formic::InputBase & userinp) {

            // reduce accumulations to root process
            { S x =  _le_avg; formic::mpi::reduce(&x,  &_le_avg, 1, MPI::SUM); }
            { S x = _le2_avg; formic::mpi::reduce(&x, &_le2_avg, 1, MPI::SUM); }
            if ( _vw > 0.0 ) {
              formic::mpi::reduce(_led_avg);
              formic::mpi::reduce(_led_le_avg);
              formic::mpi::reduce(_dr_avg);
              formic::mpi::reduce(_dr_le_avg);
              formic::mpi::reduce(_dr_le2_avg);
            }
            formic::mpi::reduce(_sref);
            formic::mpi::reduce(_href);

            // root process scales accumulations by the total weight
            if ( formic::mpi::rank() == 0 ) {
              _le_avg /= formic::unity(S()) * this->_total_weight * this->_vgsa;
              _le2_avg /= formic::unity(S()) * this->_total_weight * this->_vgsa;
              if ( _vw > 0.0 ) {
                _led_avg /= formic::unity(S()) * this->_total_weight * this->_vgsa;
                _led_le_avg /= formic::unity(S()) * this->_total_weight * this->_vgsa;
                _dr_avg /= formic::unity(S()) * this->_total_weight * this->_vgsa;
                _dr_le_avg /= formic::unity(S()) * this->_total_weight * this->_vgsa;
                _dr_le2_avg /= formic::unity(S()) * this->_total_weight * this->_vgsa;
              }
              _sref /= formic::unity(S()) * this->_total_weight * this->_vgsa;
              _href /= formic::unity(S()) * this->_total_weight * this->_vgsa;
            }

            // root process finishes building the hamiltonian matrix
            if ( _vw > 0.0 && formic::mpi::rank() == 0 ) {

              // get some constants
              const int m = _href.rows();
              const S one = formic::unity(S());

              // get the derivatives of the energy
              formic::Vector<S,int> en_der(m, formic::zero(S()));
              en_der.axpy( 2.0 * one, _dr_le_avg);
              en_der.axpy( 1.0 * one, _led_avg);
              en_der.axpy(-2.0 * _le_avg, _dr_avg);

              // add final contributions to hamiltonian matrix
              formic::xgemm('N', 'N', m, m, 1, -1.0 * _vw * one, &en_der.at(0), m, &_led_avg.at(0), 1, one, &_href(0,0), m);
              formic::xgemm('N', 'N', m, m, 1, -1.0 * _vw * one, &_led_avg.at(0), m, &en_der.at(0), 1, one, &_href(0,0), m);
              formic::xgemm('N', 'N', m, m, 1,  1.0 * _vw * one, &en_der.at(0), m, &en_der.at(0), 1, one, &_href(0,0), m);
              for (int i = 1; i < m; i++) {
                const S d = _vw * (   1.0 * _dr_le2_avg(i)
                                    + 1.0 * _led_le_avg(i)
                                    - 1.0 * _le2_avg * _dr_avg(i)
                                    - 1.0 * _le_avg * en_der(i) );
                _href(i,0) += d;
                _href(0,i) += d;
              }
              _href.axpy(_vw * ( _le2_avg - _le_avg * _le_avg ), _sref);
              formic::of << boost::format("_le_avg = %20.12f") % _le_avg << std::endl;
              formic::of << boost::format("    var = %20.12f") % ( _le2_avg - _le_avg * _le_avg ) << std::endl;

            }

          }

          ///////////////////////////////////////////////////////////////////////////////////////////////////
          /// \brief  resets the accumulator
          ///
          ///////////////////////////////////////////////////////////////////////////////////////////////////
          void child_reset(const formic::Configuration & config) {
            _href = formic::zero(S());
            _sref = formic::zero(S());
          }

      };

    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructor
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      HamOvlpBuilderLM(const formic::InputBase & userinp,
                       boost::shared_ptr<formic::fqmc::Sampler<S> > & sampler,
                       const boost::shared_ptr<formic::CoeffObj<S> > & vf,
                       const double shift)
        : _ham(1+vf->n_der(), 1+vf->n_der(), formic::zero(S())),
          _ovl(1+vf->n_der(), 1+vf->n_der(), formic::zero(S())),
          _shift( formic::zero(S()) )
      {

        // build the matrices
        {

          formic::start_timer("LM build matrices");

          // print some info
          if ( formic::mpi::rank() == 0 )
            formic::of << boost::format("explicitly building the linear method's hamiltonian and overlap matrices") << std::endl << std::endl;

          // prepare accumulator
          const double vw = userinp.get<double>("lm_var_weight");
          boost::shared_ptr<formic::Accumulator<S> > acc( new HamOvlpAcc(sampler->gf(), vf, _ham, _ovl, vw) );

          // take the sample
          sampler->simple_sample() << vf << acc;

          formic::stop_timer("LM build matrices");

        }

        // add the hamiltonian shift
        this->update_ham_shift(shift);

        //// print
        //if ( formic::mpi::rank() == 0 ) {
        //  formic::of << boost::format("_ham") << std::endl;
        //  for (int i = 0; i < _ham.rows(); i++) {
        //    for (int j = 0; j < _ham.cols(); j++)
        //      formic::of << boost::format("  %14.6e") % _ham(i,j);
        //    formic::of << std::endl;
        //  }
        //  formic::of << std::endl;
        //  formic::of << boost::format("_ovl") << std::endl;
        //  for (int i = 0; i < _ovl.rows(); i++) {
        //    for (int j = 0; j < _ovl.cols(); j++)
        //      formic::of << boost::format("  %14.6e") % _ovl(i,j);
        //    formic::of << std::endl;
        //  }
        //  formic::of << std::endl;
        //}

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  updates the diagonal shift
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update_ham_shift(const double new_shift) {
        const S old_shift = _shift;
        _shift = formic::unity(S()) * new_shift;
        for (int i = 1; i < _ham.rows(); i++)
          _ham(i,i) += _shift - old_shift;
      }

      /// \brief  returns a reference to the (shifted) hamiltonian matrix
      formic::Matrix<S,int> & ham() { return _ham; }

      /// \brief  returns a reference to the overlap matrix
      formic::Matrix<S,int> & ovl() { return _ovl; }

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a base class for objects that compute matrix multiplication in the context of the
  ///         linear method optimizer
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class MatOpLM : public formic::MatrixActor<S> {

    private:

      /// \brief  disable copy constructor
      MatOpLM(const MatOpLM & other);

      /// \brief  disable assignment operator
      MatOpLM & operator=(const MatOpLM & other);

    protected:

      /// \brief  user input options
      const formic::InputBase & _userinp;

//      /// \brief  a walker that will be used in sampling
//      formic::VMC_Walker<S> _walker;

      /// \brief  the object that controls the sampling
      boost::shared_ptr<formic::fqmc::Sampler<S> > _sampler;

      /// \brief  the coefficient object for the value function
      boost::shared_ptr<formic::CoeffObj<S> > _vf;

      /// \brief  an accumulator that will be used in sampling
      boost::shared_ptr<formic::Accumulator<S> > _acc;

      /// \brief  the Q matrix
      formic::Matrix<S,int> _qmat;

      /// \brief  number of samples
      int _nsamp;

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  a function to do the necessary sampling for the matrix multiplication
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void do_sampling() {
        _sampler->simple_sample(_nsamp) << _vf << _acc;
      }

    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructor
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      MatOpLM(const formic::InputBase & userinp,
              boost::shared_ptr<formic::fqmc::Sampler<S> > sampler,
              const boost::shared_ptr<formic::CoeffObj<S> > & vf,
              int nsamp = -1)
        : _userinp(userinp),
          _sampler(sampler),
          _vf(vf),
          _acc(),
          _nsamp(nsamp)
      {

        if ( _nsamp < 0 )
          _nsamp = _sampler->sl();

        // compute the hamiltonian and overlap matrices (this should be replaced by a more efficient sampling of just what is needed for the Q matrix)
//        formic::fqmc::HamOvlpBuilderLM<S> hob(userinp, co, init_config, lfg, nsamp, 0.0);
//      
//        // build the Q matrix, which projects the current wave function out of the wave function derivatives
//        if ( formic::mpi::rank() == 0 ) {
//          const int n = 1 + vf->n_der();
//          _qmat.reset(n, n, formic::zero(S()));
//          for (int i = 0; i < n; i++) {
//            _qmat.at(0,i) = -hob.ovl().at(0,i);
//            _qmat.at(i,i) = formic::unity(S());
//          }
//          //formic::of << _qmat.print("%14.6e", "_qmat") << std::endl;
//        }

      }

      formic::Matrix<S,int> & qmat() { return _qmat; }

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class to multiply by the hamiltonian in the context of the linear method optimizer
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class HamOpLM : public formic::fqmc::MatOpLM<S> {

    private:

      /// \brief  the shift to apply to the diagonal of the Hamiltonian
      S _shift;

    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructor
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      HamOpLM(const formic::InputBase & userinp,
              boost::shared_ptr<formic::fqmc::Sampler<S> > & sampler,
              const boost::shared_ptr<formic::CoeffObj<S> > & vf,
              const double shift,
              int nsamp = -1)
        : MatOpLM<S>(userinp, sampler, vf, nsamp),
          _shift( formic::unity(S()) * shift )
      {
        boost::shared_ptr<formic::Accumulator<S> > temp_acc( new formic::fqmc::DerivMultHamAcc<S>(userinp, sampler->gf(), vf) );
        this->_acc = temp_acc;
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  function to return the name of the child class
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      std::string class_name() { return "formic::fqmc::HamOpLM"; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  updates the diagonal shift
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update_ham_shift(const double new_shift) { _shift = formic::unity(S()) * new_shift; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  function used to multiply by the Hamiltonian  ( y = H d )
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void operate(const int n, S * const d, S * const y, const int myrank, const int tid) {

        // put a matrix around the input vector
        formic::Matrix<S,int> dlocal(n, 1, d);

        //// multiply by Q matrix
        //if ( formic::mpi::rank() == 0 )
        //  dlocal = this->_qmat * dlocal; // after this, root's dlocal is a new, separate matrix

        // send the vector to all processes
        formic::mpi::bcast((S *)&dlocal.at(0,0), n);

        // do the sampling
        formic::fqmc::DerivMultHamAcc<S> * ptr = dynamic_cast<formic::fqmc::DerivMultHamAcc<S> *>(this->_acc.get());
        if (ptr == 0)
          throw formic::Exception("dynamic cast returned zero pointer in HamOpLM::operate");
        ptr->set_dc(&dlocal.at(0,0));
        ptr->set_vp(y);
        this->do_sampling();

        //// multiply by Q matrix transpose
        //if ( formic::mpi::rank() == 0 )
        //  ( this->_qmat.t() * formic::ConstMatrix<S,int>(n, 1, y) ) >>= y;

        // add shift
        for (int i = 1; i < n; i++)
          y[i] += _shift * d[i];

        // send the vector to all processes
        formic::mpi::bcast((S *)y, n);

      }

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class to multiply by the overlap matrix in the context of the linear method optimizer
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class OvlpOpLM : public formic::fqmc::MatOpLM<S> {

    private:

      /// \brief  the shift to apply to the diagonal of the overlap matrix
      S _shift;

    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructor
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      OvlpOpLM(const formic::InputBase & userinp,
               boost::shared_ptr<formic::fqmc::Sampler<S> > sampler,
               const boost::shared_ptr<formic::CoeffObj<S> > & vf,
               int nsamp = -1)
        : MatOpLM<S>(userinp, sampler, vf, nsamp),
          _shift( formic::unity(S()) * userinp.get<double>("lm_s_diag_shift") )
      {
        boost::shared_ptr<formic::Accumulator<S> > temp_acc( new formic::fqmc::DerivMultOvlpAcc<S>(userinp, sampler->gf(), vf) );
        this->_acc = temp_acc;
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  function to return the name of the child class
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      std::string class_name() { return "formic::fqmc::OvlpOpLM"; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  function used to multiply by the overlap matrix  ( y = S d )
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void operate(const int n, S * const d, S * const y, const int myrank, const int tid) {

        // put a matrix around the input vector
        formic::Matrix<S,int> dlocal(n, 1, d);

        //// multiply by Q matrix
        //if ( formic::mpi::rank() == 0 )
        //  dlocal = this->_qmat * dlocal; // after this, root's dlocal is a new, separate matrix

        // send the vector to all processes
        formic::mpi::bcast((S *)&dlocal.at(0,0), n);

        // do the sampling
        formic::fqmc::DerivMultOvlpAcc<S> * ptr = dynamic_cast<formic::fqmc::DerivMultOvlpAcc<S> *>(this->_acc.get());
        if (ptr == 0)
          throw formic::Exception("dynamic cast returned zero pointer in OvlpOpLM::operate");
        ptr->set_dc(&dlocal.at(0,0));
        ptr->set_vp(y);
        this->do_sampling();

        //// multiply by Q matrix transpose
        //if ( formic::mpi::rank() == 0 )
        //  ( this->_qmat.t() * formic::ConstMatrix<S,int>(n, 1, y) ) >>= y;

        // add shift
        for (int i = 0; i < n; i++)
          y[i] += _shift * d[i];

        // send the vector to all processes
        formic::mpi::bcast((S *)y, n);

      }

  };

//  ///////////////////////////////////////////////////////////////////////////////////////////////////
//  /// \brief  a class to accumulate linear method preconditioner data
//  ///
//  ///////////////////////////////////////////////////////////////////////////////////////////////////
//  template <class S> class LinMethPCAcc : public formic::Accumulator<S> {
//
//    private:
//
//      /// \brief  pointer to the coefficient object to use
//      boost::shared_ptr<formic::CoeffObj<S> > _co;
//
//      /// \brief  the first row of the approximate Hamiltonian
//      std::vector<S> _hrow;
//
//      /// \brief  the diagonal of the approximate Hamiltonian
//      std::vector<S> _hdiag;
//
//      /// \brief  the first row of the overlap matrix
//      std::vector<S> _srow;
//
//      /// \brief  the diagonal of the overlap matrix
//      std::vector<S> _sdiag;
//
//      /// \brief  number of preconditioner indices for each variable
//      formic::Vector<int,int> _pc_n_inds;
//
//      /// \brief  preconditioner indices
//      formic::Vector<int,int> _pc_inds;
//
//      /// \brief  preconditioner hamiltonian values
//      formic::Vector<S,int> _pc_vals_h;
//
//      /// \brief  preconditioner overlap values
//      formic::Vector<S,int> _pc_vals_s;
//
//      /// \brief  a vector to hold derivative ratios
//      std::vector<S> _der_ratios;
//
////      /// \brief  the shift to apply to the diagonal of the Hamiltonian
////      S _h_shift;
////
////      /// \brief  the shift to apply to the diagonal of the overlap matrix
////      S _s_shift;
//
//    public:
//
//      /// \brief  return the accumulator's name
//      std::string name() const { return "LinMethPCAcc accumulator"; }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief  constructs the accumulator
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      LinMethPCAcc(const formic::InputBase & userinp,
//                   const boost::shared_ptr<formic::CoeffObj<S> > & co)
//        : _co(co),
//          _hrow (1+co->n_der(), formic::zero(S())),
//          _hdiag(1+co->n_der(), formic::zero(S())),
//          _srow (1+co->n_der(), formic::zero(S())),
//          _sdiag(1+co->n_der(), formic::zero(S())),
//          _pc_n_inds(1+co->n_der(), 0),
//          _pc_inds(co->n_pc_nonzero(), 0),
//          _pc_vals_h(co->n_pc_nonzero(), formic::zero(S())),
//          _pc_vals_s(co->n_pc_nonzero(), formic::zero(S())),
//          _der_ratios(1+co->n_der(), formic::zero(S()))
//          //_h_shift( formic::unity(S()) * h_shift ),
//          //_s_shift( formic::unity(S()) * userinp.get<double>("lm_s_diag_shift") )
//      {
//        formic::start_timer("LM prep precon sparsity");
//        co->set_pc_inds(&_pc_n_inds[0], &_pc_inds[0]);
//        formic::stop_timer("LM prep precon sparsity");
//      }
//
//      /// \brief  return the first row of the approximate Hamiltonian
//      const std::vector<S> & hrow() { return _hrow; }
//
//      /// \brief  return the diagonal of the approximate Hamiltonian
//      const std::vector<S> & hdiag() { return _hdiag; }
//
//      /// \brief  return the first row of the overlap matrix
//      const std::vector<S> & srow() { return _srow; }
//
//      /// \brief  return the diagonal of the overlap matrix
//      const std::vector<S> & sdiag() { return _sdiag; }
//
//      /// \brief  return the numbers of preconditioning indices
//      const formic::Vector<int,int> & pc_n_inds() { return _pc_n_inds; }
//
//      /// \brief  return the preconditioning indices
//      const formic::Vector<int,int> & pc_inds() { return _pc_inds; }
//
//      /// \brief  return the preconditioning h values
//      const formic::Vector<S,int> & pc_vals_h() { return _pc_vals_h; }
//
//      /// \brief  return the preconditioning h values
//      const formic::Vector<S,int> & pc_vals_s() { return _pc_vals_s; }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief Accumulate the local contribution.
//      ///
//      /// \param[in]      userinp  the user's input options
//      /// \param[in]      w        the weight to apply to this configuration's contribution
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void child_accumulate(const formic::InputBase & userinp, const double w) {
//
//        // get the local derivative ratios
//        _der_ratios.at(0) = formic::unity(S());
//        _co->der_ratios(&_der_ratios.at(0)+1);
//
//        // get the local diagonal energy
//        const S lde = _co->local_diag_energy();
//
//        // update the preconditioner values
//        {
//          const int * i_ptr = &_pc_inds[0];
//          S * h_ptr = &_pc_vals_h[0];
//          S * s_ptr = &_pc_vals_s[0];
//          for (int i = 0; i < _pc_n_inds.size(); i++) {
//            if ( _der_ratios.at(i) != formic::zero(S()) ) {
//              const int * stop_ptr = i_ptr + _pc_n_inds[i];
//              const S pre_h = w * lde * formic::conj(_der_ratios[i]);
//              const S pre_s = w * formic::conj(_der_ratios[i]);
//              for ( ; i_ptr != stop_ptr; i_ptr++, h_ptr++, s_ptr++) {
//                const S dr = _der_ratios[*i_ptr];
//                *h_ptr += pre_h * dr;
//                *s_ptr += pre_s * dr;
//              }
//            } else {
//              i_ptr += _pc_n_inds[i];
//              h_ptr += _pc_n_inds[i];
//              s_ptr += _pc_n_inds[i];
//            }
//          }
//        }
//
//        // update the first hamiltonian row
//        formic::xaxpy(_der_ratios.size(), w * lde, &_der_ratios.at(0), 1, &_hrow.at(0), 1);
//
//        // update the first overlap matrix row
//        formic::xaxpy(_der_ratios.size(), w, &_der_ratios.at(0), 1, &_srow.at(0), 1);
//
//        // get the square norms of the derivative ratios
//        for (int i = 0; i < _der_ratios.size(); i++)
//          _der_ratios[i] = formic::conj(_der_ratios[i]) * _der_ratios[i];
//
//        // update the hamiltonian diagonal
//        formic::xaxpy(_der_ratios.size(), w * lde, &_der_ratios.at(0), 1, &_hdiag.at(0), 1);
//
//        // update the overlap matrix diagonal
//        formic::xaxpy(_der_ratios.size(), w, &_der_ratios.at(0), 1, &_sdiag.at(0), 1);
//
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief finalizes the accumulation (this may include for example communicating the results)
//      ///
//      /// \param[in]      userinp  the user's input options
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void child_finalize(const formic::InputBase & userinp) {
//
//        // reduce preconditioner h array
//        if (formic::mpi::rank() == 0) {
//          formic::Vector<S,int> red_src = _pc_vals_h.clone();
//          formic::mpi::reduce(&red_src[0], &_pc_vals_h[0], _pc_inds.size(), MPI::SUM);
//          formic::xscal(_pc_inds.size(), formic::unity(S()) / this->_total_weight, &_pc_vals_h[0], 1);
//        } else {
//          formic::mpi::reduce(&_pc_vals_h[0], (S*)0, _pc_inds.size(), MPI::SUM);
//        }
//
//        // reduce preconditioner s array
//        if (formic::mpi::rank() == 0) {
//          formic::Vector<S,int> red_src = _pc_vals_s.clone();
//          formic::mpi::reduce(&red_src[0], &_pc_vals_s[0], _pc_inds.size(), MPI::SUM);
//          formic::xscal(_pc_inds.size(), formic::unity(S()) / this->_total_weight, &_pc_vals_s[0], 1);
//        } else {
//          formic::mpi::reduce(&_pc_vals_s[0], (S*)0, _pc_inds.size(), MPI::SUM);
//        }
//
//        // reduce hamiltonian row
//        formic::xcopy(_der_ratios.size(), &_hrow.at(0), 1, &_der_ratios.at(0), 1);
//        formic::mpi::reduce(&_der_ratios.at(0), &_hrow.at(0), _der_ratios.size(), MPI::SUM);
//        formic::xscal(_der_ratios.size(), formic::unity(S()) / this->_total_weight, &_hrow.at(0), 1);
//
//        // reduce hamiltonian diagonal
//        formic::xcopy(_der_ratios.size(), &_hdiag.at(0), 1, &_der_ratios.at(0), 1);
//        formic::mpi::reduce(&_der_ratios.at(0), &_hdiag.at(0), _der_ratios.size(), MPI::SUM);
//        formic::xscal(_der_ratios.size(), formic::unity(S()) / this->_total_weight, &_hdiag.at(0), 1);
//
//        // reduce overlap row
//        formic::xcopy(_der_ratios.size(), &_srow.at(0), 1, &_der_ratios.at(0), 1);
//        formic::mpi::reduce(&_der_ratios.at(0), &_srow.at(0), _der_ratios.size(), MPI::SUM);
//        formic::xscal(_der_ratios.size(), formic::unity(S()) / this->_total_weight, &_srow.at(0), 1);
//
//        // reduce overlap diagonal
//        formic::xcopy(_der_ratios.size(), &_sdiag.at(0), 1, &_der_ratios.at(0), 1);
//        formic::mpi::reduce(&_der_ratios.at(0), &_sdiag.at(0), _der_ratios.size(), MPI::SUM);
//        formic::xscal(_der_ratios.size(), formic::unity(S()) / this->_total_weight, &_sdiag.at(0), 1);
//
//        //// add shifts
//        //for (int i = 1; i < _hdiag.size(); i++)
//        //  _hdiag[i] += _h_shift;
//        //for (int i = 0; i < _sdiag.size(); i++)
//        //  _sdiag[i] += _s_shift;
//
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief  resets the accumulator
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void child_reset(const formic::Configuration & config) {
//        std::fill(_hrow.begin(), _hrow.end(), formic::zero(S()));
//        std::fill(_hdiag.begin(), _hdiag.end(), formic::zero(S()));
//        std::fill(_srow.begin(), _srow.end(), formic::zero(S()));
//        std::fill(_sdiag.begin(), _sdiag.end(), formic::zero(S()));
//        _pc_vals_h = formic::zero(S());
//        _pc_vals_s = formic::zero(S());
//      }
//
//  };

//  ///////////////////////////////////////////////////////////////////////////////////////////////////
//  /// \brief  a class used to operate by the linear method preconditioner
//  ///
//  ///////////////////////////////////////////////////////////////////////////////////////////////////
//  template <class S> class PrecOpLM : public formic::MatrixActor<S> {
//
//    private:
//
//      /// \brief  accumulator used to compute the preconditioner
//      boost::shared_ptr<formic::Accumulator<S> > _acc;
//
//      /// \brief  derived class pointer to the accumulator
//      formic::fqmc::LinMethPCAcc<S> * _acc_ptr;
//
//      /// \brief  a reference to the current approximation to the eigenvalue in the davidson solver
//      const double & _eval;
//
//      /// \brief  max number of generalized minimum residual iterations
//      int _gmr_iter;
//
//      /// \brief  number generalized minimum residual iterations before restarting
//      int _gmr_rest;
//
//      /// \brief  the shift to apply to the diagonal of the Hamiltonian
//      S _h_shift;
//
//      /// \brief  the shift to apply to the diagonal of the overlap matrix
//      S _s_shift;
//
//      /// \brief  work space
//      std::vector<S> _work;
//
//    private:
//
//      void form_row_and_diagonal(S * const row, S * const diag) const {
//
//        // get pointers to the hamiltonian row and diagonal and the overlap row and diagonal
//        const S * const hr = &_acc_ptr->hrow().at(0);
//        const S * const hd = &_acc_ptr->hdiag().at(0);
//        const S * const sr = &_acc_ptr->srow().at(0);
//        const S * const sd = &_acc_ptr->sdiag().at(0);
//
//        const int n = _acc_ptr->hrow().size();
//
//        // form row
//        for (int i = 0; i < n; i++)
//          row[i] = hr[i] - _eval * sr[i];
//        row[0] -= _eval * _s_shift;
//
//        // form diagonal
//        const double zero_thresh = 1.0e-11;
//        const double small_thresh = 1.0e-8;
//        S h_shift_to_use = formic::zero(S());
//        for (int i = 0; i < n; i++) {
//          const S val = ( hd[i] + h_shift_to_use ) - _eval * ( sd[i] + _s_shift );
//          const double abs_val = std::abs(abs_val);
//          if      ( abs_val < zero_thresh  )
//            diag[i] = formic::unity(S()) * small_thresh;
//          else if ( abs_val < small_thresh )
//            diag[i] = small_thresh * ( val / abs_val );
//          else
//            diag[i] = val;
//          h_shift_to_use = _h_shift;
//        }
//
//      }
//
//    public:
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief  constructor
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      PrecOpLM(const formic::InputBase & userinp,
//               const boost::shared_ptr<formic::CoeffObj<S> > & co,
//               const formic::Configuration & init_config,
//               const formic::LaggedFibonacci<> & lfg,
//               const double & eval,
//               const double h_shift)
//        : _eval(eval),
//          _gmr_iter(userinp.get<int>("lm_genminres_iter")),
//          _gmr_rest(userinp.get<int>("lm_genminres_rest")),
//          _acc_ptr(0),
//          _h_shift( formic::unity(S()) * h_shift ),
//          _s_shift( formic::unity(S()) * userinp.get<double>("lm_s_diag_shift") ),
//          _work(1+co->n_der(), formic::zero(S()))
//      {
//
//        // prepare the accumulator
//        {
//          typedef boost::shared_ptr<formic::Accumulator<S> > ShAcPtr;
//          ShAcPtr temp_acc( new formic::fqmc::LinMethPCAcc<S>(userinp, co) );
//          _acc = temp_acc;
//        }
//
//        // get a derived class pointer to the accumulator
//        _acc_ptr = dynamic_cast<formic::fqmc::LinMethPCAcc<S> *>(_acc.get());
//        if (_acc_ptr == 0)
//          throw formic::Exception("dynamic cast returned zero pointer in PrecOpLM::PrecOpLM");
//
//        // prepare a walker
//        formic::VMC_Walker<S> walker(userinp);
//        walker.gf() = co;
//        walker.accumulators().push_back(_acc);
//        walker.config() = init_config;
//        walker.lfg() = lfg;
//
//        // accumulate hamiltonian and overlap matrix row and column
//        formic::start_timer("LM precon elements");
//        if ( !userinp.get<bool>("exact_sampling") )
//          formic::sample(userinp, userinp.get<int>("opt_samp_len"), walker);
//        else
//          walker.exact_sample(userinp);
//        formic::stop_timer("LM precon elements");
//
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief  function to return the name of the child class
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      std::string class_name() { return "formic::fqmc::PrecOpLM"; }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief  updates the diagonal shift
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void update_ham_shift(const double new_shift) { _h_shift = formic::unity(S()) * new_shift; }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief  operates by the matrix A in A x = b, where x is the result of preconditioning on b
//      ///
//      /// \param[in]      v      the vector on which A will operate
//      ///
//      /// \return  the matrix-vector product A * v
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      formic::Vector<S,int> operator*(const formic::Vector<S,int> & v) const {
//
//        // initialize vector to return
//        formic::Vector<S,int> retval(v.size(), formic::zero(S()));
//
//        // row by row, multiply by the matrix A
//        const int * i_ptr = &_acc_ptr->pc_inds()[0];
//        const S * h_ptr = &_acc_ptr->pc_vals_h()[0];
//        const S * s_ptr = &_acc_ptr->pc_vals_s()[0];
//        const formic::Vector<int,int> & pc_n_inds = _acc_ptr->pc_n_inds();
//        for (int i = 0; i < pc_n_inds.size(); i++) {
//          assert( pc_n_inds[i] >= 0 );
//          const int * stop_ptr = i_ptr + pc_n_inds[i];
//          for ( ; i_ptr != stop_ptr; i_ptr++, h_ptr++, s_ptr++)
//            retval[i] += ( (*h_ptr) - _eval * (*s_ptr) ) * v[*i_ptr];
//        }
//
//        // add the shifts
//        retval[0] -= _eval * _s_shift * v[0];
//        for (int i = 1; i < retval.size(); i++)
//          retval[i] += ( _h_shift - _eval * _s_shift ) * v[i];
//
//        // return the result
//        return retval;
//
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief  solves the pre-preconditioning equation used in genminres
//      ///
//      /// \param[in]      v      the vector to pre-precondition
//      ///
//      /// \return  the pre-preconditioned vector
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      formic::Vector<S,int> solve(const formic::Vector<S,int> & v) const {
//
//        // get and check the dimension
//        const int n = _acc_ptr->hdiag().size();
//        if ( v.size() != n )
//          throw formic::Exception("dimension mismatch in formic::fqmc::PrecOpLM::solve");
//
//        // initialize vector to return
//        formic::Vector<S,int> retval = v.clone();
//
//        // re-form the combined row and combined diagonal
//        formic::Vector<S,int> row(n, formic::zero(S()));
//        formic::Vector<S,int> diag(n, formic::zero(S()));
//        this->form_row_and_diagonal(&row[0], &diag[0]);
//
//        // gaussian elimination of first row
//        S upL = diag[0];
//        for (int i = 1; i < n; i++) {
//          upL -= formic::square_norm(row[i]) / diag[i];
//          retval[0] -= retval[i] * row[i] / diag[i];
//        }
//
//        // avoid dividing by zero
//        if ( std::abs(upL) < 1.0e-11 )
//          upL = formic::unity(S()) * 1.0e-6;
//        else if ( std::abs(upL) < 1.0e-6 )
//          upL = 1.0e-6 * ( upL / std::abs(upL) );
//
//        // gaussian elimination of first column
//        for (int i = 1; i < n; i++)
//          retval[i] -= retval[0] * formic::conj(row[i]) / upL;
//
//        // diagonal division
//        retval[0] = retval[0] / upL;
//        for (int i = 1; i < n; i++)
//          retval[i] = retval[i] / diag[i];
//
//        // return the result
//        return retval;
//
//      }
//
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      /// \brief  function used to multiply by the inverse of the approximate matrix pencil
//      ///
//      /// y = ( H_approx - eval * S_approx )^(-1) d
//      ///
//      ///////////////////////////////////////////////////////////////////////////////////////////////////
//      void operate(const int n, S * const d, S * const y, const int myrank, const int tid) {
//
//        formic::start_timer("LM preconditioning");
//
//        if ( formic::mpi::rank() == 0 ) {
//
//          int m = _gmr_rest;
//          int mi = _gmr_iter;
//          S tol = formic::unity(S()) * 1.0e-8;
//          formic::Matrix<S,int> temp_mat(m+1, m+1, formic::zero(S()));
//          formic::genminres_classes::VectorNorm<S,int> norm;
//          formic::genminres_classes::VectorDot<S,int> dot;
//          formic::Vector<S,int> x(n, formic::zero(S()));
//          formic::Vector<S,int> b(n, formic::zero(S()));
//          formic::xcopy(n, d, 1, &b[0], 1);
//          formic::genminres(*this,
//                            x,
//                            b,
//                            *this,
//                            temp_mat,
//                            norm,
//                            dot,
//                            m,
//                            mi,
//                            tol);
//          formic::xcopy(n, &x[0], 1, y, 1);
//
//        }
//
//        formic::stop_timer("LM preconditioning");
//
//      }
//
//  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class used to operate by a short-sample linear method preconditioner
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class ShortSampPrecOpLM : public formic::MatrixActor<S> {

    private:

      /// \brief  a reference to the current approximation to the eigenvalue in the davidson solver
      const double & _eval;

      /// \brief  max number of generalized minimum residual iterations
      int _gmr_iter;

      /// \brief  number generalized minimum residual iterations before restarting
      int _gmr_rest;

      /// \brief  the shift to apply to the diagonal of the Hamiltonian
      S _h_shift;

      /// \brief  the shift to apply to the diagonal of the overlap matrix
      S _s_shift;

      /// \brief  operator that performs multiplication by the Hamiltonian
      formic::fqmc::HamOpLM<S> _opH;

      /// \brief  operator that performs multiplication by the overlap matrix
      formic::fqmc::OvlpOpLM<S> _opS;

      /// \brief  work space
      std::vector<S> _work;

    private:

    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructor
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      ShortSampPrecOpLM(const formic::InputBase & userinp,
                        boost::shared_ptr<formic::fqmc::Sampler<S> > & sampler,
                        const boost::shared_ptr<formic::CoeffObj<S> > & vf,
                        const int nsamp,
                        const double & eval,
                        const double h_shift)
        : _eval(eval),
          _gmr_iter(userinp.get<int>("lm_genminres_iter")),
          _gmr_rest(userinp.get<int>("lm_genminres_rest")),
          _h_shift( formic::unity(S()) * h_shift ),
          _s_shift( formic::unity(S()) * userinp.get<double>("lm_s_diag_shift") ),
          _opH(userinp, sampler, vf, h_shift, nsamp),
          _opS(userinp, sampler, vf, nsamp),
          _work(1+vf->n_der(), formic::zero(S()))
      {}

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  function to return the name of the child class
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      std::string class_name() { return "formic::fqmc::ShortSampPrecOpLM"; }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  updates the diagonal shift
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update_ham_shift(const double new_shift) {
        _h_shift = formic::unity(S()) * new_shift;
        _opH.update_ham_shift(new_shift);
      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  operates by the matrix A in A x = b, where x is the result of preconditioning on b
      ///
      /// \param[in]      v      the vector on which A will operate
      ///
      /// \return  the matrix-vector product A * v
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      formic::Vector<S,int> operator*(const formic::Vector<S,int> & v) {

        // initialize vector to return
        formic::Vector<S,int> retval(v.size(), formic::zero(S()));

        _opH.operate(_work.size(), &v.clone()[0], &retval[0], formic::mpi::rank(), 0);
        _opS.operate(_work.size(), &v.clone()[0], &_work[0], formic::mpi::rank(), 0);
        formic::xaxpy(_work.size(), -_eval * formic::unity(S()), &_work[0], 1, &retval[0], 1);
        formic::mpi::bcast(&retval[0], retval.size());

        // return the result
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  solves the pre-preconditioning equation used in genminres
      ///
      /// \param[in]      v      the vector to pre-precondition
      ///
      /// \return  the pre-preconditioned vector
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      formic::Vector<S,int> solve(const formic::Vector<S,int> & v) const {

        // return the result (currently no pre-preconditioner is being used)
        return v.clone();

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  function used to multiply by the inverse of the approximate matrix pencil
      ///
      /// y = ( H_approx - eval * S_approx )^(-1) d
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void operate(const int n, S * const d, S * const y, const int myrank, const int tid) {

        formic::start_timer("LM preconditioning");

        formic::mpi::bcast((S *)(d), n);

        if ( true || formic::mpi::rank() == 0 ) {

          int m = _gmr_rest;
          int mi = _gmr_iter;
          S tol = formic::unity(S()) * 1.0e-8;
          formic::Matrix<S,int> temp_mat(m+1, m+1, formic::zero(S()));
          formic::genminres_classes::VectorNorm<S,int> norm;
          formic::genminres_classes::VectorDot<S,int> dot;
          formic::Vector<S,int> x(n, formic::zero(S()));
          formic::Vector<S,int> b(n, formic::zero(S()));
          formic::xcopy(n, d, 1, &b[0], 1);
          formic::genminres(*this,
                            x,
                            b,
                            *this,
                            temp_mat,
                            norm,
                            dot,
                            m,
                            mi,
                            tol);
          formic::xcopy(n, &x[0], 1, y, 1);

        }

        formic::stop_timer("LM preconditioning");

      }

  };

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief  a class used as the davidson solver for the linear method
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class S> class DavidsonLM {

    // data members
    private:

      /// \brief  number of wave function derivatives
      int _nder;

      /// \brief  dimension of the first derivative subspace
      int _nfds;

      /// \brief  dimension of the krylov subspace
      int _nkry;

      /// \brief  maximum number of iterations to use when solving the eigenproblem
      int _max_iter;

      /// \brief  residual threshold to use when solving the eigenproblem
      double _thresh;

      /// \brief  the current estimate of the Hamiltonian eigenvalue
      double _energy;

      /// \brief  energy weighting
      double _ew;

      /// \brief  variance weighting
      double _vw;

      /// \brief  initial wavefunction energy estimate
      double _init_energy;

      /// \brief  initial wavefunction standard deviation estimate
      double _init_sdev;

      /// \brief  maximum change allowed from initial energy
      double _max_energy_change;

      /// \brief  shift applied to the Hamiltonian diagonal
      double _hshift;

      /// \brief  the smallest subspace overlap eigenvalue
      double _min_S_eval;

      /// \brief  residual threshold above which to use arnoldi updates
      double _arnoldi_thresh;

      /// \brief  flag to tell whether the lowest eigenvalue has a complex component
      bool _eval_was_complex;

      /// \brief  operator that performs multiplication by the Hamiltonian
      boost::shared_ptr<formic::fqmc::HamOpLM<S> > _opH;

      /// \brief  operator that performs multiplication by the overlap matrix
      boost::shared_ptr<formic::fqmc::OvlpOpLM<S> > _opS;

      /// \brief  operator that performs multiplication by the preconditioner
      boost::shared_ptr<formic::fqmc::ShortSampPrecOpLM<S> > _opM;
      //boost::shared_ptr<formic::fqmc::PrecOpLM<S> > _opM;

      /// \brief  class that builds the hamiltonian and overlap matrices
      boost::shared_ptr<formic::fqmc::HamOvlpBuilderLM<S> > _hs_builder;

      /// \brief  the (normalized) krylov subspace basis vectors
      std::vector<std::vector<S> > _kvecs;

      /// \brief  vectors resulting from Hamiltonian action on the krylov basis vectors
      std::vector<std::vector<S> > _hvecs;

      /// \brief  vectors resulting from overlap action on the krylov basis vectors
      std::vector<std::vector<S> > _svecs;

      /// \brief  subspace projection of the Hamiltonian
      std::vector<S> _subH;

      /// \brief  subspace projection of the overlap
      std::vector<S> _subS;

      /// \brief  work vector
      std::vector<S> _wv1;

      /// \brief  work vector
      std::vector<S> _wv2;

      /// \brief  subspace eigenvector
      std::vector<S> _sub_evec;

      /// \brief  overlap matrix eigenvectors used in direct diagonalization
      formic::Matrix<S,int> _ovl_evecs;

      /// \brief  overlap matrix eigenvalues used in direct diagonalization
      formic::Vector<double,int> _ovl_evals;

    // private functions
    private:

      class NBiggerThan {
        private:
          double _x;
        public:
          NBiggerThan(double x) : _x(x) {}
          int operator()(int total, double val) const { return total + ( val > _x ? 1 : 0 ); }
      };

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  solves the subspace eigenvalue problem
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void solve_subspace() {

        const int m = _nkry;

        // create arrays and other variables used in diagonalization routines
        std::vector<S> matA(m*m, formic::zero(S()));
        std::vector<S> matB(m*m, formic::zero(S()));
        std::vector<double> wr(m, 0.0);
        std::vector<double> wi(m, 0.0);
        std::vector<S> wc(m, formic::zero(S()));
        int lwork = (10 + m) * m;
        std::vector<S> work(lwork, formic::zero(S()));
        std::vector<double> rwork(lwork, 0.0);
        int info = 0;

        // diagonalize the subspace S matrix
        std::vector<S> s_evecs(_subS);
        //// !!! BEGIN BUG CHECKING !!!
        //formic::of << formic::Matrix<S,int>(m, m, &s_evecs.at(0)).print("%10.2e", "krylov subspace overlap matrix") << std::endl;
        //// !!!  END BUG CHECKING  !!!
        std::vector<double> s_evals(m, 0.0);
        formic::xsyev('V', 'U', m, &s_evecs.at(0), m, &s_evals.at(0),
                      &work.at(0), lwork, &rwork.at(0), info);
        if ( info != 0 )
          throw formic::Exception("formic::xsyev failed for subspace S matrix with info = %i in formic::DavidsonLM::solve_subspace") % info;

        // get the smallest eigenvalue of the subspace S matrix
        _min_S_eval = *std::min_element(s_evals.begin(), s_evals.end());

        // compute the product subS^(-1) subH
        formic::xgemm('T', 'N', m, m, m, formic::unity(S()), &s_evecs.at(0), m, &_subH.at(0), m, formic::zero(S()), &matB.at(0), m);
        for (int j = 0; j < m; j++)
        for (int i = 0; i < m; i++)
          matB.at(i+m*j) /= s_evals.at(i);
        formic::xgemm('N', 'N', m, m, m, formic::unity(S()), &s_evecs.at(0), m, &matB.at(0), m, formic::zero(S()), &matA.at(0), m);

        // compute the eigenvalues and right eigenvectors of subS^(-1) subH
        formic::xgeev('N', 'V', m, &matA.at(0), m, &wc.at(0), &wr.at(0), &wi.at(0), &matB.at(0), m, &matB.at(0), m,
                      &work.at(0), lwork, &rwork.at(0), info);
        if ( info != 0 )
          throw formic::Exception("formic::xgeev failed for  subS^(-1) subH  with info = %i in formic::DavidsonLM::solve_subspace") % info;

        // find the lowest eigenpair
        double min_wr = wr.at(0);
        double min_wi = wi.at(0);
        int min_i = 0;
        for (int i = 0; i < m; i++) {
          if ( wr.at(i) < min_wr ) {
            min_wr = wr.at(i);
            min_wi = wi.at(i);
            min_i = i;
          }
        }

        // set the energy
        _energy = min_wr;

        // if the lowest eigenvalue has an imaginary component, we abort
        _eval_was_complex = false;
        if ( std::abs(min_wi) > 1.0e-10 ) {
          _eval_was_complex = true;
          return;
        }

        // if the lowest eigenvalue is real, use it and its corresponding eigenvector
        _sub_evec.resize(m);
        for (int i = 0; i < m; i++)
          _sub_evec.at(i) = matB.at( i + m * min_i );

        // normalize the eigenvector
        S evec_norm = formic::zero(S());
        for (int i = 0; i < m; i++)
        for (int j = 0; j < m; j++)
          evec_norm += formic::conj(_sub_evec.at(i)) * _sub_evec.at(j) * _subS.at(i+m*j);
        formic::xscal(m, formic::unity(S()) / std::abs(std::sqrt(evec_norm)), &_sub_evec.at(0), 1);

      }

    // public functions
    public:

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  constructor
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      DavidsonLM(const formic::InputBase & userinp,
                 const boost::shared_ptr<formic::CoeffObj<S> > & vf,
                 const double init_energy,
                 const double init_sdev)
        : _nder(vf->n_der()),
          _nfds(1+vf->n_der()),
          _nkry(0),
          _max_iter(userinp.get<int>("lm_davidson_iter")),
          _thresh(userinp.get<double>("lm_davidson_thresh")),
          _energy(0.0),
          _ew(1.0 - userinp.get<double>("lm_var_weight")),
          _vw(userinp.get<double>("lm_var_weight")),
          _init_energy(init_energy),
          _init_sdev(init_sdev),
          _max_energy_change(userinp.get<double>("lm_max_e_change")),
          _hshift(0.0),
          _min_S_eval(1.0),
          _arnoldi_thresh(userinp.get<double>("lm_arnoldi_thresh")),
          _wv1(1+vf->n_der(), formic::zero(S())),
          _wv2(1+vf->n_der(), formic::zero(S()))
      {}

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  resets the solver
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void reset(const formic::InputBase & userinp,
                 boost::shared_ptr<formic::fqmc::Sampler<S> > & sampler,
                 const boost::shared_ptr<formic::CoeffObj<S> > & vf,
                 const double shift,
                 const bool rebuild_matrices = false)
      {

        // save the shift
        _hshift = shift;

        // empty the krylov subspace
        _kvecs.clear();
        _hvecs.clear();
        _svecs.clear();
        _subH.clear();
        _subS.clear();
        _sub_evec.clear();
        _min_S_eval = 1.0;
        _nkry = 0;

        // empty the overlap eigenvectors and eigenvalues
        _ovl_evecs.clear();
        _ovl_evals.clear();

        // initializations for matrix-multiplication-by-sampling
        if ( ! userinp.get<bool>("lm_build_matrices") ) {

          // initialize the Hamiltonian actor
          _opH.reset();
          boost::shared_ptr<formic::fqmc::HamOpLM<S> > tempH( new formic::fqmc::HamOpLM<S>(userinp, sampler, vf, shift) );
          _opH = tempH;

          // initialize the overlap actor
          _opS.reset();
          boost::shared_ptr<formic::fqmc::OvlpOpLM<S> > tempS( new formic::fqmc::OvlpOpLM<S>(userinp, sampler, vf) );
          _opS = tempS;

          // initialize the preconditioner actor
          _opM.reset();
          boost::shared_ptr<formic::fqmc::ShortSampPrecOpLM<S> > tempM( new formic::fqmc::ShortSampPrecOpLM<S>(userinp, sampler, vf, userinp.get<int>("pc_samp_len"), _energy, shift) );
          //boost::shared_ptr<formic::fqmc::PrecOpLM<S> > tempM( new formic::fqmc::PrecOpLM<S>(userinp, gf, vf, init_config, lfg, _energy, shift) );
          _opM = tempM;

        // initializations for building matrices directly
        } else {

          // if necessary, initialize the hamiltonian and overlap matrix builder (this computes the matrices by sampling)
          if (rebuild_matrices)
            _hs_builder.reset();
          if ( ! _hs_builder ) {
            boost::shared_ptr<formic::fqmc::HamOvlpBuilderLM<S> > tempHSB( new formic::fqmc::HamOvlpBuilderLM<S>(userinp, sampler, vf, shift) );
            _hs_builder = tempHSB;
          }

        }

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  adds a new krylov basis vector
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void add_krylov_vector(const formic::InputBase & userinp, const formic::Vector<S,int> & v) {

        // check vector length
        if ( formic::mpi::rank() == 0 && v.size() != _nfds )
          throw formic::Exception("bad vector length of %i in DavidsonLM::add_krylov_vector:  expected length of %i") % v.size() % _nfds;

        //// prepare work vectors
        //_wv1.resize(_nfds);
        //_wv2.resize(_nfds);
        //std::fill(_wv1.begin(), _wv1.end(), formic::zero(S()));
        //std::fill(_wv2.begin(), _wv2.end(), formic::zero(S()));

        // increment krylov subspace size and remember old size
        const int nold = _nkry++;

        // multiply new vector by the overlap
        if (formic::mpi::rank() == 0)
          formic::xcopy(_nfds, &v.at(0), 1, &_wv1.at(0), 1);
        if ( ! userinp.get<bool>("lm_build_matrices") ) {
          _opS->operate(_nfds, &_wv1.at(0), &_wv2.at(0), formic::mpi::rank(), 0);
        } else if (formic::mpi::rank() == 0) {
          formic::Matrix<S,int>(_nfds, 1, &_wv2.at(0)) <<= _hs_builder->ovl() * formic::Matrix<S,int>(_nfds, 1, &_wv1.at(0));
        }

        //// normalize the new vector and save it and the overlap's action upon in
        //if (formic::mpi::rank() == 0) {
        //  const S norm_inv = formic::unity(S()) / std::abs(std::sqrt(formic::xdotc(_nfds, &_wv1.at(0), 1, &_wv2.at(0), 1)));
        //  formic::xscal(_nfds, norm_inv, &_wv1.at(0), 1);
        //  formic::xscal(_nfds, norm_inv, &_wv2.at(0), 1);
        //  _kvecs.push_back(_wv1);
        //  _svecs.push_back(_wv2);
        //}

        // normalize and graham schmidt the new vector and save it and the overlap's action upon in
        if (formic::mpi::rank() == 0) {
          //// !!! BEGIN BUG CHECKING !!!
          //formic::of << boost::format("new krylov vector norm before graham schmidt = %16.4e")
          //                % std::abs(std::sqrt(formic::xdotc(_nfds, &_wv1.at(0), 1, &_wv2.at(0), 1)))
          //           << std::endl;
          //// !!!  END BUG CHECKING  !!!
          const S norm_inv = formic::unity(S()) / std::abs(std::sqrt(formic::xdotc(_nfds, &_wv1.at(0), 1, &_wv2.at(0), 1)));
          formic::xscal(_nfds, norm_inv, &_wv1.at(0), 1);
          for (int i = 0; i < nold; i++)
            formic::xaxpy(_nfds, -formic::xdotc(_nfds, &_svecs.at(i).at(0), 1, &_wv1.at(0), 1), &_kvecs.at(i).at(0), 1, &_wv1.at(0), 1);
          //// !!! BEGIN BUG CHECKING !!!
          //formic::of << formic::Matrix<S,int>(_nfds, 1, &_wv1.at(0)).print("%10.2e", "update vector after orthogonalization");
          //// !!!  END BUG CHECKING  !!!
        }
        if ( ! userinp.get<bool>("lm_build_matrices") ) {
          _opS->operate(_nfds, &_wv1.at(0), &_wv2.at(0), formic::mpi::rank(), 0);
        } else if (formic::mpi::rank() == 0) {
          formic::Matrix<S,int>(_nfds, 1, &_wv2.at(0)) <<= _hs_builder->ovl() * formic::Matrix<S,int>(_nfds, 1, &_wv1.at(0));
        }
        if (formic::mpi::rank() == 0) {
          //// !!! BEGIN BUG CHECKING !!!
          //formic::of << formic::Matrix<S,int>(_nfds, 1, &_wv2.at(0)).print("%10.2e", "S * vector");
          //const S xdotc_result = formic::xdotc(_nfds, &_wv1.at(0), 1, &_wv2.at(0), 1);
          //formic::of << boost::format("xdotc result = %16.4e") % xdotc_result << std::endl;
          //if ( formic::real(xdotc_result) < 0.0 ) {
          //  const int m = _nfds;
          //  //std::vector<S> matA(m*m, formic::zero(S()));
          //  //std::vector<S> matB(m*m, formic::zero(S()));
          //  //std::vector<double> wr(m, 0.0);
          //  //std::vector<double> wi(m, 0.0);
          //  //std::vector<S> wc(m, formic::zero(S()));
          //  formic::Matrix<double,int> w(m, 1, 0.0);
          //  int lwork = (10 + m) * m;
          //  std::vector<S> work(lwork, formic::zero(S()));
          //  std::vector<double> rwork(lwork, 0.0);
          //  int info = 0;
          //  formic::Matrix<S,int> s_evecs = _hs_builder->ovl().clone();
          //  formic::xsyev('V', 'U', m, &s_evecs.at(0,0), m, &w.at(0,0), &work.at(0), lwork, &rwork.at(0), info);
          //  if ( info != 0 )
          //    throw formic::Exception("formic::xsyev failed for S matrix with info = %i in my check for positive definite overlap") % info;
          //  formic::of << w.print("%10.2e", "overall S eigenvalues");
          //}
          //formic::of << boost::format("new krylov vector norm  after graham schmidt = %16.4e")
          //                % std::abs(std::sqrt(formic::xdotc(_nfds, &_wv1.at(0), 1, &_wv2.at(0), 1)))
          //           << std::endl;
          //// !!!  END BUG CHECKING  !!!
          const S norm_inv = formic::unity(S()) / std::abs(std::sqrt(formic::xdotc(_nfds, &_wv1.at(0), 1, &_wv2.at(0), 1)));
          formic::xscal(_nfds, norm_inv, &_wv1.at(0), 1);
          formic::xscal(_nfds, norm_inv, &_wv2.at(0), 1);
          _kvecs.push_back(_wv1);
          _svecs.push_back(_wv2);
        }

        // multiply new vector by the Hamiltonian
        if ( ! userinp.get<bool>("lm_build_matrices") ) {
          _opH->operate(_nfds, &_wv1.at(0), &_wv2.at(0), formic::mpi::rank(), 0);
        } else if (formic::mpi::rank() == 0) {
          //_hs_builder->update_ham_shift(_hshift);
          formic::Matrix<S,int>(_nfds, 1, &_wv2.at(0)) <<= _hs_builder->ham() * formic::Matrix<S,int>(_nfds, 1, &_wv1.at(0));
        }

        if (formic::mpi::rank() == 0) {

          // save new Hamiltonian action vector
          _hvecs.push_back(_wv2);

          // update subspace projection of the Hamiltonian
          std::vector<S> newH(_nkry*_nkry, formic::zero(S()));
          for (int i = 0; i < nold; i++)
          for (int j = 0; j < nold; j++)
            newH.at(i+_nkry*j) = _subH.at(i+nold*j);
          _subH = newH;
          _subH.at(_nkry-1+_nkry*(_nkry-1)) = formic::xdotc(_nfds, &_kvecs.at(_nkry-1).at(0), 1, &_hvecs.at(_nkry-1).at(0), 1);
          for (int i = 0; i < nold; i++) {
            _subH.at(_nkry-1+_nkry*i        ) = formic::xdotc(_nfds, &_kvecs.at(_nkry-1).at(0), 1, &_hvecs.at(i      ).at(0), 1);
            _subH.at(i      +_nkry*(_nkry-1)) = formic::xdotc(_nfds, &_kvecs.at(i      ).at(0), 1, &_hvecs.at(_nkry-1).at(0), 1);
          }

          // update subspace projection of the overlap
          std::vector<S> newS(_nkry*_nkry, formic::zero(S()));
          for (int i = 0; i < nold; i++)
          for (int j = 0; j < nold; j++)
            newS.at(i+_nkry*j) = _subS.at(i+nold*j);
          _subS = newS;
          _subS.at(_nkry-1+_nkry*(_nkry-1)) = formic::xdotc(_nfds, &_kvecs.at(_nkry-1).at(0), 1, &_svecs.at(_nkry-1).at(0), 1);
          for (int i = 0; i < nold; i++) {
            _subS.at(_nkry-1+_nkry*i        ) = formic::xdotc(_nfds, &_kvecs.at(_nkry-1).at(0), 1, &_svecs.at(i      ).at(0), 1);
            _subS.at(i      +_nkry*(_nkry-1)) = formic::xdotc(_nfds, &_kvecs.at(i      ).at(0), 1, &_svecs.at(_nkry-1).at(0), 1);
          }

        }

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  changes the Hamiltonian shift
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      void update_ham_shift(const formic::InputBase & userinp, const double new_shift) {

        // shift update for solver with matrix multiplication by sampling
        if ( ! userinp.get<bool>("lm_build_matrices") ) {

          _opH->update_ham_shift(new_shift);
          _opM->update_ham_shift(new_shift);

        // shift update for direct matrix building
        } else {

          _hs_builder->update_ham_shift(new_shift);

        }

        // shift update for krylov vectors
        if ( formic::mpi::rank() == 0 && new_shift != _hshift ) {

          for (int i = 0; i < _nkry; i++)
            formic::xaxpy(_nder, formic::unity(S()) * ( new_shift - _hshift ), &_kvecs.at(i).at(1), 1, &_hvecs.at(i).at(1), 1);

          for (int j = 0; j < _nkry; j++)
          for (int i = 0; i < _nkry; i++)
            _subH.at(i+_nkry*j) = formic::xdotc(_nfds, &_kvecs.at(i).at(0), 1, &_hvecs.at(j).at(0), 1);

        }

        // remember the shift
        _hshift = new_shift;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  solves the eigenvalue problem via the davidson method
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      bool iterative_solve(const formic::InputBase & userinp, double & eval, formic::Vector<S,int> & evec) {

        // return value:  whether the davidson/arnoldi iterations satisfied the various perscribed stability criteria
        bool retval = true;

        // iterate
        bool converged = false;
        int iter = 0;
        while (true) {

          // solve subspace eigenvalue problem
          if (formic::mpi::rank() == 0)
            this->solve_subspace();
          formic::mpi::bcast(_energy);
          eval = _energy;

          // check if energy had an imaginary component
          formic::mpi::bcast(_eval_was_complex);
          if ( _eval_was_complex ) {
            retval = false;
            if (formic::mpi::rank() == 0)
              formic::of << boost::format("davidson iteration %4i:  stopping due to imaginary component in energy") % iter << std::endl;
            break;
          }

          // check if energy change is still reasonable
          if ( std::abs( _energy - _ew * _init_energy - _vw * _init_sdev * _init_sdev ) > _max_energy_change ) {
            retval = false;
            if (formic::mpi::rank() == 0)
              formic::of << boost::format("davidson iteration %4i:  stopping due to too-large eigenvalue change") % iter << std::endl;
            break;
          }

          // if the overlap matrix has become singular, stop
          formic::mpi::bcast(_min_S_eval);
          if ( _min_S_eval < userinp.get<double>("lm_min_S_eval") ) {
            if (formic::mpi::rank() == 0) {
              formic::of << boost::format("davidson iteration %4i:  stopping due to small subspace S eigenvalue of %.2e")
                            % iter
                            % _min_S_eval
                         << std::endl;
              //throw formic::Exception("AAAAAAAHHHH  small subspace overlap eigenvalue OOOHHHH NOOOOOOOO");
            }
            break;
          }

          // construct the new eigenvector estimate
          if (formic::mpi::rank() == 0) {
            evec = formic::zero(S());
            for (int i = 0; i < _nkry; i++)
              formic::xaxpy(_nfds, _sub_evec.at(i), &_kvecs.at(i).at(0), 1, &evec.at(0), 1);
          }

          // construct the new residual
          //_wv1.resize(_nfds);
          std::fill(_wv1.begin(), _wv1.end(), formic::zero(S()));
          if (formic::mpi::rank() == 0) {
            for (int i = 0; i < _nkry; i++) {
              formic::xaxpy(_nfds, _sub_evec.at(i), &_hvecs.at(i).at(0), 1, &_wv1.at(0), 1);
              formic::xaxpy(_nfds, -_energy * _sub_evec.at(i), &_svecs.at(i).at(0), 1, &_wv1.at(0), 1);
            }
          }

          //if (formic::mpi::rank() == 0) {
          //  formic::of << std::endl;
          //  formic::of << boost::format("printing residual:") << std::endl;
          //  for (int i = 0; i < _wv1.size(); i++)
          //    formic::of << boost::format("%20.12f") % _wv1.at(i) << std::endl;
          //  formic::of << std::endl;
          //}

          // Here we compute the norm of the residual <Psi^i| ( H - E ) |Psi> / <Psi|Psi>.
          // Note that the "normalization", i.e. dividing by <Psi|Psi>, is done automatically by the
          // markov chain monte carlo sampling, so we don't need to multiply the residual by the
          // overlap matrix in order to compute it's norm.
          double rms_res = 1.0e100;
          //_opS->operate(_nfds, &_wv1.at(0), &_wv2.at(0), formic::mpi::rank(), 0);
          //if (formic::mpi::rank() == 0)
          //  rms_res = std::abs(std::sqrt( formic::xdotc(_nfds, &_wv1.at(0), 1, &_wv2.at(0), 1) / double(_nfds) ));
          if (formic::mpi::rank() == 0)
            rms_res = std::abs(std::sqrt( formic::xdotc(_nfds, &_wv1.at(0), 1, &_wv1.at(0), 1) / double(_nfds) ));
          formic::mpi::bcast(rms_res);

          // print iteration result
          if (formic::mpi::rank() == 0)
            formic::of << boost::format("davidson iteration %4i:    krylov dim = %3i   eval = %20.12f        rms_res = %.2e      min_S_eval = %.2e")
                          % iter
                          % _nkry
                          % _energy
                          % rms_res
                          % _min_S_eval
                          //% ( need_to_symmetrize ? "    symm" : "" )
                          //% ( E > last_eval      ? "    up" : "" )
                       << std::endl;

          // check for convergence
          converged = ( rms_res < _thresh );
          if (converged)
            break;

          // increment the iteration count and stop if the maximum number of iterations has been reached
          if ( ++iter >= _max_iter )
            break;

          // compute the new subspace basis vector
          if ( rms_res > _arnoldi_thresh ) { // arnoldi when residual is large
            if (formic::mpi::rank() == 0)
              formic::xcopy(_nfds, &_wv1.at(0), 1, &_wv2.at(0), 1);
          } else { // davidson when residual is small
            if ( userinp.get<bool>("lm_build_matrices") )
              throw formic::Exception("no preconditioner implemented for explicitly built matrices in formic::DavidsonLM::iterative_solve");
            _opM->operate(_nfds, &_wv1.at(0), &_wv2.at(0), formic::mpi::rank(), 0);
          }

          // expand the krylov subspace
          formic::Vector<S,int> to_add;
          to_add <<= _wv2;
          this->add_krylov_vector(userinp, to_add);

        }

        // print iteration result
        if (converged && formic::mpi::rank() == 0)
          formic::of << boost::format("davidson solver converged in %10i iterations") % iter << std::endl << std::endl;
        else if (formic::mpi::rank() == 0)
          formic::of << boost::format("davidson solver did not converge after %10i iterations") % iter << std::endl << std::endl;

        //// multiply by the Q matrix
        //if (formic::mpi::rank() == 0)
        //  ( _opS->qmat() * formic::ConstMatrix<S,int>(evec.size(), 1, &evec.at(0)) ) >>= evec;

        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  solves the eigenvalue problem
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      bool solve(const formic::InputBase & userinp, double & eval, formic::Vector<S,int> & evec) {

        // check evec length
        if ( evec.size() != _nfds )
          throw formic::Exception("input eigenvector has size %i (should be %i) in formic::DavidsonLM::solve") % evec.size() % _nfds;

        // if requested, solve iteratively via the davidson method
        //if ( userinp.get<bool>("lm_do_davidson") )
        if ( true )
          return this->iterative_solve(userinp, eval, evec);

        formic::start_timer("LM direct diagonalization");

        // return value:  whether the iterations stayed within the perscribed energy limits
        bool retval = true;

        // root process performs diagonalization
        if ( formic::mpi::rank() == 0 ) {

          //////////////////////////////////////////
          //   adapted from solve subspace code   //
          //////////////////////////////////////////

          // get dimension
          const int n = _hs_builder->ham().rows();

          // create arrays and other variables used in diagonalization routines
          formic::Vector<double,int> wr(n, 0.0);
          formic::Vector<double,int> wi(n, 0.0);
          formic::Vector<S,int> wc(n, formic::zero(S()));
          int lwork = (10 + n) * n;
          formic::Vector<S,int> work(lwork, formic::zero(S()));
          formic::Vector<double,int> rwork(lwork, 0.0);
          int info = 0;

//          // diagonalize the overlap matrix and store its evals and evecs if we have not done so already
//          if ( _ovl_evecs.rows() != n || _ovl_evecs.cols() != n ) {
//
//            _ovl_evecs <<= _hs_builder->ovl();
//            _ovl_evals.reset(n, 0.0); // = formic::Vector<double,int>(n, 0.0);
//            formic::xsyev('V', 'U', n, &_ovl_evecs(0,0), n, &_ovl_evals.at(0), &work.at(0), lwork, &rwork.at(0), info);
//            if ( info != 0 )
//              throw formic::Exception("formic::xsyev failed for overlap matrix with info = %i in formic::DavidsonLM::solve") % info;
//
//          }
//
//          // check that the overlap eigenvalues are positive
//          if ( _ovl_evals.at(n-1) <= 0.0 )
//            throw formic::Exception("overlap matrix eigenvalues are all negative in formic::DavidsonLM::solve");
//
//          // print info on S eigenvalues
//          formic::of << boost::format("min overlap eigenvalue = %20.6e") % _ovl_evals.at(0)   << std::endl;
//          formic::of << boost::format("max overlap eigenvalue = %20.6e") % _ovl_evals.at(n-1) << std::endl;
//          formic::of << std::endl;

          // for convenience, construct a matrix around the output eigenvector
          formic::Matrix<S,int> full_evec(n, 1, &evec.at(0));

          // build Q matrix, which projects the current wave function out of the wave function derivatives
          formic::Matrix<S,int> qmat(n, n, formic::zero(S()));
          for (int i = 0; i < n; i++) {
            qmat.at(0,i) = -_hs_builder->ovl().at(0,i);
            qmat.at(i,i) = formic::unity(S());
          }

          //// make sure the shift is up to date
          //_hs_builder->update_ham_shift(_hshift);

          //// print hamiltonian matrix
          //formic::of << ( _hs_builder->ham() - _hs_builder->ham().c() ).print("%14.6e", "ham deviation from symmetric") << std::endl;

          //// print overlap matrix
          //formic::of << ( _hs_builder->ovl() - _hs_builder->ovl().c() ).print("%14.6e", "ovl deviation from symmetric") << std::endl;

          //// print Q matrix
          //formic::of << qmat.print("%14.6e", "qmat") << std::endl;

          // get Q-modified hamiltonian and overlap matrices
          formic::Matrix<S,int> qham = qmat.c() * _hs_builder->ham() * qmat;
          formic::Matrix<S,int> qovl = qmat.c() * _hs_builder->ovl() * qmat; // + 1.0e-13 * formic::unity(S()) * qmat.c() * qmat;

          //// print Q-modified overlap matrix
          //formic::of << qovl.print("%14.6e", "qovl") << std::endl;

          // diagonalize Q-modified overlap matrix
          formic::Vector<double,int> qovl_evals(n, 0.0);
          formic::xsyev('V', 'U', n, &qovl.at(0,0), n, &qovl_evals.at(0), &work.at(0), lwork, &rwork.at(0), info);
          if ( info != 0 )
            throw formic::Exception("formic::xsyev failed for Q-modified overlap matrix with info = %i in formic::DavidsonLM::solve") % info;

          // print info on overlap eigenvalues
          formic::of << boost::format("min overlap eigenvalue = %13.6e") % qovl_evals.at(0)   << std::endl;
          formic::of << boost::format("max overlap eigenvalue = %13.6e") % qovl_evals.at(n-1) << std::endl;
          formic::of << std::endl;

          // get overlap cutoff
          //const S overlap_shift = formic::unity(S()) * (-2.0);
          //const S overlap_shift = formic::unity(S()) * userinp.get<double>("lm_s_diag_shift");
          const double overlap_shift = userinp.get<double>("lm_s_diag_shift") * qovl_evals.at(n-1);
          formic::of << boost::format("overlap eigenvectors with eigenvalues below %.6e will be truncated") % overlap_shift << std::endl << std::endl;

          // count the number of eigenvalues below the threshold
          const int p = std::accumulate(qovl_evals.begin(), qovl_evals.end(), int(0), NBiggerThan(overlap_shift));

          //for (int i = 0; i < n; i++)
          //  formic::of << boost::format("qovl eval %4i  =  %14.6e") % i % qovl_evals.at(i) << std::endl;
          //formic::of << std::endl;

          //// print Q-modified overlap matrix eigenvectors
          //formic::of << qovl.print("%14.6e", "qovl") << std::endl;

          //formic::of << boost::format("discarded %i basis functions") % (n-p) << std::endl;
          //formic::of << std::endl;

          // for each eigenvalue above the cutoff, get its square root and its corresponding eigenvector
          formic::Vector<double,int> trunc_sqrts(p, 0.0);
          formic::Matrix<S,int> trunc_evecs(n, p, formic::zero(S()));
          for (int j = n-p, q = 0; j < n; j++, q++) {
            trunc_sqrts.at(q) = std::abs(std::sqrt(qovl_evals.at(j)));
            for (int i = 0; i < n; i++)
              trunc_evecs.at(i,q) = qovl.at(i,j);
          }

          // get truncated Hamiltonian
          formic::Matrix<S,int> trunc_ham = trunc_evecs.c() * qham * trunc_evecs;

          // scale by eigenvalue square roots
          for (int j = 0; j < p; j++)
          for (int i = 0; i < p; i++)
            trunc_ham.at(i,j) /= trunc_sqrts.at(i) * trunc_sqrts.at(j);

          //formic::of << boost::format("original ham is %i by %i:") % _hs_builder->ham().rows() % _hs_builder->ham().cols() << std::endl;
          //for (int i = 0; i < n; i++) {
          //  for (int j = 0; j < n; j++)
          //    formic::of << boost::format("  %14.6e") % _hs_builder->ham().at(i,j);
          //  formic::of << std::endl;
          //}
          //formic::of << std::endl;
          //formic::of << std::endl;

          //formic::of << boost::format("trunc_ham is %i by %i:") % trunc_ham.rows() % trunc_ham.cols() << std::endl;
          //for (int i = 0; i < p; i++) {
          //  for (int j = 0; j < p; j++)
          //    formic::of << boost::format("  %14.6e") % trunc_ham.at(i,j);
          //  formic::of << std::endl;
          //}
          //formic::of << std::endl;
          //formic::of << std::endl;

          //// save a copy of the truncated hamiltonian
          //formic::Matrix<S,int> trunc_ham_copy = trunc_ham.clone();

          //// print truncated hamiltonian
          //formic::of << trunc_ham.print("%14.6e", "trunc_ham") << std::endl;

          // compute the eigenvalues and right eigenvectors of the truncated hamiltonian
          formic::Matrix<S,int> matB(p, p, formic::zero(S()));
          formic::xgeev('N', 'V', p, &trunc_ham.at(0,0), p, &wc.at(0), &wr.at(0), &wi.at(0), &matB.at(0,0), p, &matB.at(0,0), p,
                        &work.at(0), lwork, &rwork.at(0), info);
          if ( info != 0 )
            throw formic::Exception("formic::xgeev failed for truncated hamiltonian with info = %i in formic::DavidsonLM::solve") % info;

          // find the lowest eigenpair and make sure its eigenvalue is real
          double min_wr = wr.at(0);
          double min_wi = wi.at(0);
          int min_i = 0;
          for (int i = 0; i < p; i++) {
            if ( wr.at(i) < min_wr ) {
              min_wr = wr.at(i);
              min_wi = wi.at(i);
              min_i = i;
            }
          }

          // if the lowest eigenvalue is real, use it and its corresponding eigenvector
          if ( std::abs(min_wi) < 1.0e-10 ) {
            eval = min_wr;
          } else {
            throw formic::Exception("lowest eigenvalue had an imaginary component of %.2e in formic::DavidsonLM::solve") % min_wi;
          }

//            // diagonalize truncated Hamiltonian
//            formic::Vector<double,int> h_eval(p, 0.0);
//            formic::xsyev('V', 'U', p, &trunc_ham.at(0,0), p, &h_eval.at(0), &work.at(0), lwork, &rwork.at(0), info);
//            if ( info != 0 )
//              throw formic::Exception("formic::xsyev failed for truncated hamiltonian matrix with info = %i in formic::DavidsonLM::solve") % info;
//
//            // get eigenvalue
//            eval = h_eval.at(0);

          //formic::Matrix<S,int> inner_resid_vec = trunc_ham_copy * matB - eval * formic::unity(S()) * matB;
          //double inner_rms_res = 0.0;
          //for (int i = 0; i < p; i++)
          //  inner_rms_res += formic::square_norm( inner_resid_vec(i,min_i) );
          //inner_rms_res = std::abs( std::sqrt( inner_rms_res / double(_nfds) ) );
          //formic::of << boost::format("inner_rms_res = %14.6e") % inner_rms_res << std::endl << std::endl;

          // get the lowest eigenvector divided by the square roots of the overlap eigenvalues
          formic::Matrix<S,int> scaled_h_evec(p, 1, formic::zero(S()));
          for (int i = 0; i < p; i++)
            scaled_h_evec.at(i,0) = matB.at(i,min_i) / trunc_sqrts.at(i);

          // get eigenvector in un-truncated space
          full_evec <<= qmat * trunc_evecs * scaled_h_evec;
          //full_evec <<= trunc_evecs * scaled_h_evec;

          // get the residual
          formic::Matrix<S,int> res_vec = _hs_builder->ham() * full_evec - eval * formic::unity(S()) * ( _hs_builder->ovl() * full_evec );

          // print the residual norm
          formic::of << boost::format("rms_res = %14.6e") % std::abs( std::sqrt( ( res_vec.c() * res_vec ).at(0,0) / double(_nfds) ) ) << std::endl << std::endl;
          //formic::of << res_vec.print("%14.6e", "residual") << std::endl;

//            // scale by overlap eigenvalues
//            for (int j = 0; j < p; j++)
//            for (int i = 0; i < p; i++)
//              trunc_ham.at(i,j) /= trunc_sqrts.at(i) * trunc_sqrts.at(i);
//
//            // compute the eigenvalues and right eigenvectors
//            formic::Matrix<S,int> matB(p, p, formic::zero(S()));
//            formic::xgeev('N', 'V', p, &trunc_ham.at(0,0), p, &wc.at(0), &wr.at(0), &wi.at(0), &matB(0,0), p, &matB(0,0), p,
//                          &work.at(0), lwork, &rwork.at(0), info);
//            if ( info != 0 )
//              throw formic::Exception("formic::xgeev failed for trunc_ham with info = %i in formic::DavidsonLM::solve") % info;
//  
//            // find the lowest eigenpair and make sure its eigenvalue is real
//            double min_wr = wr.at(0);
//            double min_wi = wi.at(0);
//            int min_i = 0;
//            for (int i = 0; i < p; i++) {
//              if ( wr.at(i) < min_wr ) {
//                min_wr = wr.at(i);
//                min_wi = wi.at(i);
//                min_i = i;
//              }
//            }
//  
//            // if the lowest eigenvalue is real, use it and its corresponding eigenvector
//            if ( std::abs(min_wi) < 1.0e-10 ) {
//              eval = min_wr;
//              formic::Matrix<S,int> h_evec(p, 1, formic::zero(S()));
//              for (int i = 0; i < p; i++)
//                h_evec(i,0) = matB(i,min_i);
//              formic::Matrix<S,int> full_evec(n, 1, &evec.at(0));
//              full_evec <<= trunc_evecs * h_evec;
//            } else {
//              throw formic::Exception("lowest eigenvalue had an imaginary component of %.2e in formic::DavidsonLM::solve") % min_wi;
//            }
//  
//            // normalize the eigenvector
//            S evec_norm = formic::zero(S());
//            for (int i = 0; i < n; i++)
//            for (int j = 0; j < n; j++)
//              evec_norm += formic::conj(evec[i]) * evec[j] * _hs_builder->ovl().at(i,j);
//            formic::xscal(n, formic::unity(S()) / std::abs(std::sqrt(evec_norm)), &evec.at(0), 1);

//          if (true) {
//
//          } else {
//
//            // get overlap shift
//            const S overlap_shift = formic::unity(S()) * std::max( 1.0e-9, userinp.get<double>("lm_s_diag_shift") * _ovl_evals.at(n-1) );
//            formic::of << boost::format("overlap eigenvalues will be shifted by %.4e") % formic::real(overlap_shift) << std::endl << std::endl;
//  
//            //// get the smallest eigenvalue of the subspace S matrix
//            //min_S_eval = *std::min_element(wr.begin(), wr.end());
//  
//            // compute the product S^(-1) H and store it in matA
//            formic::Matrix<S,int> matA = _hs_builder->ham().clone(); // deep copy
//            formic::xgemm('C', 'N', n, n, n, formic::unity(S()), &_ovl_evecs(0,0), n, &matA(0,0), n, formic::zero(S()), &work.at(0), n);
//            for (int j = 0; j < n; j++)
//            for (int i = 0; i < n; i++)
//              work(i+n*j) /= _ovl_evals(i) + overlap_shift;
//            formic::xgemm('N', 'N', n, n, n, formic::unity(S()), &_ovl_evecs(0,0), n, &work.at(0), n, formic::zero(S()), &matA(0,0), n);
//  
//            // compute the eigenvalues and right eigenvectors of S^(-1) H and store the evecs in matB
//            formic::Matrix<S,int> matB(n, n, formic::zero(S()));
//            formic::xgeev('N', 'V', n, &matA(0,0), n, &wc.at(0), &wr.at(0), &wi.at(0), &matB(0,0), n, &matB(0,0), n,
//                          &work.at(0), lwork, &rwork.at(0), info);
//            if ( info != 0 )
//              throw formic::Exception("formic::xgeev failed for  S^(-1) H  with info = %i in formic::DavidsonLM::solve") % info;
//  
//            // find the lowest eigenpair and make sure its eigenvalue is real
//            double min_wr = wr.at(0);
//            double min_wi = wi.at(0);
//            int min_i = 0;
//            for (int i = 0; i < n; i++) {
//              if ( wr.at(i) < min_wr ) {
//                min_wr = wr.at(i);
//                min_wi = wi.at(i);
//                min_i = i;
//              }
//            }
//  
//            // if the lowest eigenvalue is real, use it and its corresponding eigenvector
//            if ( std::abs(min_wi) < 1.0e-10 ) {
//              eval = min_wr;
//              for (int i = 0; i < n; i++)
//                evec[i] = matB(i,min_i);
//            } else {
//              throw formic::Exception("lowest eigenvalue had an imaginary component of %.2e in formic::DavidsonLM::solve") % min_wi;
//            }
//  
//          }

          // normalize the eigenvector
          S evec_norm = formic::zero(S());
          for (int i = 0; i < n; i++)
          for (int j = 0; j < n; j++)
            evec_norm += formic::conj(evec[i]) * evec[j] * _hs_builder->ovl().at(i,j);
          formic::xscal(n, formic::unity(S()) / std::abs(std::sqrt(evec_norm)), &evec.at(0), 1);

          //// print eigenvector
          //formic::of << full_evec.print("%30.12e", "evec elements") << std::endl;

          // print the result
          formic::of << boost::format("linear method direct diagonalization gave an eigenvalue of %20.12f") % eval << std::endl << std::endl;
          //formic::of << boost::format("linear method direct diagonalization gave an eigenvalue of %20.12f for the %.4e shift") % eval % _hshift << std::endl << std::endl;

        }

        // broadcast final eigenvalue
        formic::mpi::bcast(eval);

        formic::stop_timer("LM direct diagonalization");

        // return whether or not the solve was successful
        return retval;

      }

      ///////////////////////////////////////////////////////////////////////////////////////////////////
      /// \brief  returns the minimum overlap eigenvalue in the krylov subspace
      ///
      ///////////////////////////////////////////////////////////////////////////////////////////////////
      double min_S_eval() const { return _min_S_eval; }

  };

} // end namespace fqmc

} // end namespace formic

#endif
